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ABSTRACT 

The mass of the dark matter halo of the Milky Way can be estimated by fitting analytical mod¬ 
els to the phase-space distribution of dynamical tracers. We test this approach using realistic 
mock stellar halos constructed from the Aquarius N-body simulations of dark matter halos in 
the ACDM cosmology. We extend the standard treatment to include a Navarro-Frenk-White 
(NFW) potential and use a maximum likelihood method to recover the parameters describing 
the simulated halos from the positions and velocities of their mock halo stars. We find that 
the estimate of halo mass is highly correlated with the estimate of halo concentration. The 
best-fit halo masses within the virial radius, R 200 , are biased, ranging from a 40% under¬ 
estimate to a 5% overestimate in the best case (when the tangential velocities of the tracers 
are included). There are several sources of bias. Deviations from dynamical equilibrium can 
potentially cause significant bias; deviations from spherical symmetry are relatively less im¬ 
portant. Fits to stars at different galactocentric radii can give different mass estimates. By 
contrast, the model gives good constraints on the mass within the half-mass radius of tracers 
even when restricted to tracers within 60 kpc. The recovered velocity anisotropies of tracers, 
/?, are biased systematically, but this does not affect other parameters if tangential velocity 
data are used as constraints. 
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1 INTRODUCTION 

Our Milky Way (MW) galaxy provides a wealth of information on 
the physics of galaxy formation and the nature of the dark mat¬ 
ter. This information can, in principle, be unlocked from studies of 
the positions, velocities and chemistry of stars in the Galaxy, its 
satellites and globular clusters, which can be observed with high 
precision. 


bedded in it, and there are many different ways of constraining the 
MW dark matter halo masFl- 

These met h ods include timing argument estimators 
dKahn & Woltj er) 11959)) calibrated against TV-body simula- 
tions jLi_&\Vhite 2008): modeling of local cosmic expansion 


1 Penarrubia et al. 12014); the kinematics of bright satellite s 


2007bllat Barber etalj 2014 : Cautun etalj 120141) 


1 Sales et all _ _ _ _, _ __ 

particul arly Leo I jBovlan-Kolchin et alj|201 3h an d the Magellanic 


Many inferences derived from the properties of the Milky 

Clouds (Busha et al. 

2011; Gonzalez et al. 2013!); the kinematics 

Way (MW) depend on the precision and accuracy with which the 

of stellar streams 

Newbers et al. 20 id: Kiinner et al! 2015), 

mass of its dark matter halo can be estimated. An example is the 

especially the Sagittarius stream (Lawetal. 2005; Gibbons et al.. 


MW satellite galaxies with central densities as high as those of 
the most massive dark matter subhal os predicted by ACDM sim¬ 
ulations of ‘Milky Way mass’ hosts (Bovlan-Kolchi n et al JI201 lL 
l2012l:lFerrero et alj2012l ). In these simulations the number of mas¬ 
sive subhalos depends strongly on the assumed MW halo mass and 
the problem disa ppears if the MW halo mass is suffic iently small 
(<1 x 10 12 Mr Awane et all2012l:ICautun et al .1120 1 4) . 

Gravitational lensing is the most powerful method to de¬ 
termine the underlying dark matter distri butio n for lar ge sam¬ 
ples of distant galaxie s (e.g. Bartelmann & Schneider 2001; 
Mandelbaum et all 120061 : iLi e talj 120091: [Hilbert & White! [2010; 


Han et al.ll2015t) . Our MW is, however, special because we ; 


veloc i ty stars, such as those from the RAVE survey l lSmith et al] 
120071 ; IPiffl et all 120141) ; and combinations of photometric and 
kinema tic data such as Maser observations and terminal velocity 
curves JlVfc Mi 11 ar3l20 1 lLlblesti & Saluccill2013l) . Using high reso¬ 
lution hydrodynamical simulation s and the line-o f-sigh t velocity 
dispersion of tracers in the MW, iRashkov et al.l d2013l) found a 


1 We use M 200 and R 200 to denote the mass and radius of a spherical 
region with mean density equal to 200 times the critical density of the Uni¬ 
verse. 
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2 Wang et al. 


heavy MW halo mass reported in some previous measurements of 
M 200 «2x 1O 12 M 0 is disfavoured. 

Some authors have used large composite samples of objects 
assumed to be dynamical tracers in the halo, such as stars, globu¬ 
lar clusters and planetary nebulae. For example, the halo circular 
velocity, I4irc, may be inferred from the radial velocity dispersion 
of tracers, oy(r), using the spherical Jeans equation. Such meth¬ 
ods require the tracer velocity anisotropy and density profiles to 
be known or assumed. iBattaelia et al.1 <120051) made us e of a few 
hundre d stars and globular clusters from 20 to 120 kpc: IXue et al.1 
J2008h used 2401 BHB stars from SDSS/DR6 ranging from 20 to 
60 kpc; iGnedin et alJ d2010h used BHB and RR Lyrae stars rang¬ 
ing from 25 to 80 kpc; and I Watkins et al.l J2010h used 26 satellites 


within 300 k pc with tracer mass esti mators, with the method further 
improved b ylEvans et alJ J201 lh and lAn et alJd2012h . Most recently 
iKafle et all d20 1 2ll20 1 4) used a few thousand BHB stars extending 
to 60 kpc and K-giants beyond 100 kpc. 

Most measurements based on dynamical tracers involve 
assumptions about the tracer dens ity p r ofile s and velocity 
anisotropies. However, IWilkinson & Evansl d 19991) introduced a 
Bayesian likelihood analysis, based on fitting a model phase-space 
distribution function to the observed distances and velocities of 
tracers. In their analysis the tracer density profile and velocity 
anisotropy can be considered as free parameters of the distribu¬ 
tion function, to be constrained together with parameters of the 
host halo such as it s mass and characteristic sc alelength. The sam¬ 
ple of stars used bv lWilkinson & Evansl dl999i) was small and their 
best-fit host halo mass for a t runcated flat rotation_curve model was 
1.9 ^'7 x 10 12 M 0 (see also lSakamoto et al.l2003l) . More recently, 
iDeason et alJ 1 20121) used a few thousand BHB stars from SDSS up 
to r ~50 kpc. lEadie et al.l d2015h introduced a generalised Bayesian 
approach to deal with incomplete data, which avoids rewriting the 
distribution function when tangential velocities are not available. 

Fig.Usummarises the results of these studies. The x-axis is the 
measured MW halo mass. We ha ve converted results to A /200 by as¬ 
suming an NFW density profile iNavarro et ak 1996al. 1997al) and 
using the mean halo concentration relation of lDuffv et al.l d2008h in 
cases where a value for concentration is not given in the original 
study. The measurements are grouped by methodology, indicated 
by colours and labeled along the j/-axis. We group those methods 
that use large samples of dynamical tracers into two sets: 1) those 
based on the radial velocity dispersion of the tracers and spherical 
Jeans equation to infer the circular velocity and underlying poten¬ 
tial; 2) those based on fitting to model distribution functions, which 
attempt to constrain both halo mass and the velocity anisotropy of 
the tracers simultaneously. Errorbars correspond to those quoted 
by the original authors; we have converted 90% or 95% confi¬ 
dence intervals to lcr errors assuming a Gaussian distribution, ex¬ 
cept for Watkins et al. and measuremen t s in t h e PF classifica- 
tion. This is because Wilkinson & Evansl ill 9991) . Sakamoto et al.l 


120031) . Watkins et al 


| |201(]| ). and IPeason et al.1 ' 2012} ) included 
other sources of model uncertainties beyond pure statistical errors 
in their measured masses, which makes theirerrors relatively large. 
For the generalised Bayesian approach of iEadie et alJ d2015i) . we 
quote the 95% Bayesian confidence interval. Fig.[7]shows that ex¬ 
isting measurements of the most likely MW halo mass differ by 
more than a factor of 2.5, even when similar methods are used, 
although apart from a few outliers, the estimates are statistically 
consistent. 

Here we are particularly interested in methods such as that of 
IWilkinson & Evan! dl999h . which treat the spatial and dynamical 
properties of tracers as free parameters to be constrained under the 
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Figure 1 . Measured MW halo masses in the literature (ir-axis), converted 
to M 200 and categorised by methodology (y-axis). Measurements using 
similar methods and/or tracer populations are plotted with the same colour. 
Categories include the timing argument estimator (black); a model of the 
local cosmic expansion (light grey); constraints from luminous MW satel¬ 
lites such as the Magellanic Clouds (magenta), Leo I(cyan), the orbits or 
radial velocity dispersion of other bright satellites (yellow) and their Vmax 
distribution (light green); modelling of tidal streams (grey); high velocity 
stars from the RAVE survey (blue); combinations of maser observations 
and and the terminal velocity curve (pink); and (of most relevance to this 
work) dynamical modelling using large samples of dynamical tracers (red 
and green). Methods involving large samples of dynamical tracers are split 
into two categories, 1) those based on the radial velocity dispersion of trac¬ 
ers (green) and 2) those using model distribution function to constrain both 
halo properties and velocity anisotropies of tracers simultaneously (red). We 
have converted results to M 200 by assuming an NFW density profile and 
a common mass-concentration relation. 95% or 90% confidence intervals 
have been converted to lcr errors by assuming a Gaussian error distribution, 
except for Watkins et al. and measurements in the DF classification. 


assumption of theoretical phase-space distribution functions. The 
primary aim of this paper is to test the model distribution functions 
us ed in this approach. We exte nd the distribution function proposed 
by [Wilkins on & Evan s d 19991) to one based on the NFW potential, 
and model the radial profiles of tracers with a more general double 
power-law functional form. The model function is then fit to the 
phase-space distribution of stars in realistic mock stellar halo cata¬ 
logues constructed from t he cosmological gal actic halo simulations 
of the Aquarius project dSpringel et alj[200^) , to understand its re¬ 
liability and possible violations to the underlying assumptions. Our 
results have implications that are not limited to the specific form 
of the distribution function that we test, but are applicable to the 
method itself. 

This paper is structured as follows. The mock stellar halo 
catalogues are introduced in Section [2] Detailed descriptions of 
the model distribution function and the maximum likelihood ap¬ 
proach are provided in Section [3] Our results are presented in Sec- 
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tion [4] with detailed discussions of reliability and systematics in 
Section [5] and Section [6] We conclude in Section [7] Throughout 
this paper we adopt the cosmology of the Aquarius simulation se¬ 
ries (Ho = 73 kms" 1 Mpc -1 , = 0.25, Q.a = 0.75 and 
n = 1). 


2 MOCK STELLAR HALO CATALOGUE 


We use mock stellar halo catalo g ues c onstructed from the Aquarius 
N-body simulation suite (Springel et al. 2008) with the particle tag¬ 
ging method described by Cooper et al J fcoioh . to which we refer 
the reader for further details. In this section we summarise the most 
important features of these catalogues. 


2.1 The Aquarius simulations 

The Aquarius halos come from dark matter N-body simulations in 
a standard ACDM cosmology. Cos mological paramete rs are those 
from the first year data of WMAP dSpergel et al.ll2003l ). Our work 
uses the second highest resolution level of the Aquarius suite, 
which corresponds to a particle mass of ~ 10 4 /i _1 Mq. 

The simulation suite includes six dark matter halos with virial 
masses spanning the factor-of-two range of Milky Way observa¬ 
tions discussed in the previous section. We have only used five out 
of the six halos for our analysis (labeled halo A to halo E accord¬ 
ing to the Aquarius convention). The halo we have not used (halo 
F) undergoes two major merger events at z < 0.6, and is thus an 
unlikely host for a MW-like disc galaxy. We list in Table Q]the host 
halo mass, M 200 , and oth er pro perties of the five halos, which are 
taken from lNavarro et all d2010l) . 


2.2 The galaxy formation and evolution model 

The Durham semi-analytical galaxy formation model, GALFORM, 
has been used to post-process the Aquarius simulations, predicting 
the evolution of galaxies embedded in dark matter halos. To con¬ 
struct the mock stellar h alo catalo gues used in this paper, the ver¬ 
sion described by iFont et aD J20 1 II) was ado pted. This model ha s 
several mino r differences from the model of iBower et al] J2006l ). 
such that the lFont et aDHimb model matches better the observed 
luminosity function, luminosity-metallicity relation and radial dis¬ 
tribution of MW satellites. The main changes are a more self- 
consistent calculation of the effects of the patronisation background 
and a higher chemical yield in supemovae feedback. 


2.3 Particle tagging 

The GALFORM model predicts the amount of stellar mass present 
in each dark matter halo in the simulation at each output time, as 
well as properties of stellar populations such as their total metal- 
licity. However, GALFORM does not provide detailed informa¬ 
tion about how these stars are distributed in galaxies. The particle 
tagging method of ICooper et alj d2010h is a way to determine the 
six-dimensional spatial and velocity distribution of stars from dark 
matter only simulations, by associating newly-formed stars with 
tightly bound dark matter particles. 

At each simulation snapshot, each newly formed stellar popu¬ 
lation predicted by GALFORM is assigned to the 1% most bound 
dark matter particles in its host dark matter halo. Each “tagged” 
dark matter particle then represents a fraction of a single stellar 
population, the age and metallicity of which are also known from 



Figure 2. The radial density profiles of stellar mass (red points and lines 
with errors) in the mock stellar halo catalogue of the five Aquarius halos. 
The errors are in most cases of comparable size to the symbols and are 
almost invisible. Dashed black curves are double power-law fits to the red 
data points obtained from a x 2 minimisation. Green dashed curves are the 
best-fit density profiles from the maximum likelihood method. 


GALFORM. Traced forward to the present day, these tagged par¬ 
ticles give predictions for the observed luminosity functions and 
structural properties of MW and M31 satellites that match well to 
observations. Recently. iCooper et al .1 j20 l~3t) have applied this tech¬ 
nique to large-scale cosmological simulations and have shown that 
it produces galactic surface brightness profiles that agree well with 
the outer regions of stacked galaxy profiles from SDSS. 

Our study is based on tagged dark matter particles from ac¬ 
creted satellite galaxies. We ignore particles associated with in situ 
star formation in the central galaxy. Strictly, our results thus only 
apply in the case where most MW halo stars originate from ac¬ 
cretion. This is supported by the data of iBell et al] <20081 l2010h 
although other work suggests that a certain fraction of the halo 
stars are contributed by in-situ star for mation, especially c l ose to 
the central galaxy (r < 30 kpc) (s ee, e.g. lCarollo et all2007ll201o| : 
IZolotov et alj|20icl ; iHelmi et~aPl201 lh . Ignoring the possible in- 
situ component is thus a weakness of our mock stellar halo cata¬ 
logue. Nevertheless, our mock halo stars enable us to test and con¬ 
strain the theoretical distribution function and, in practice, most of 
our conclusions (see Sec.[4]and Sec. [5} do not depend on whether 
the MW halo stars formed in-situ or were brought in by accretion. 


3 METHODOLOGY 

In this section we discuss the theoretical context of our method for 
constraining dark matter halo properties using dynamical tracers 
and a maximum likelihood approach based on theoretical distribu¬ 
tion functions. In Sec. 13.11 we describe how the phase-space dis¬ 
tribution of the tracer population is modeled. Sec. l3.2| gives details 
about the explicit form of the distribution function. The likelihood 
function is introduced in Sec. 13.31 Finally, we describe how we 
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4 Wang et al. 


weight tagged particles and how errors are estimated i n Sec. 13.41 
Our method follows that of 1 Wilkinson & Evansl d 19991 ) but intro¬ 
duces significant modifications to the form of the dark matter halo 
potential and the assumed tracer density profile. 


3.1 Phase-space distribution of Milky Way halo stars 

The phase-space distribution function of tracers (e.g. stars) bound 
to a dark matter halo potential (binding energy E > 0) can be de¬ 
scribed by the Eddington formula dEddingtonll 19161) . The simplest 
isotropic and spherically symmetric case is 


FIE), 


;/ 


d p(r) d$(r) 


\/8tt 2 d E V maX|t) d$(r) y/E-^r)' 


( 1 ) 


where the distribution function only depends on the binding energy 
per unit mass, E = <f>(r) — \. <f>(r) and v 2 /2 are the underlying 
dark matter halo potential and kinetic energy per unit mass of trac¬ 
ers. The integral goes from the potential at the tracer boundar}@ to 
the binding energy of interest. Usually both the zero point of po¬ 
tential and tracer boundary, r maXj t, are chosen at infinity, and thus 
4?(r max ,t) = 0. 

In reality the velocity distribution of tracers may be 
anisotropic and depend both on energy and angular momentum, 
L. In the simplest case, the distribution function is assumed to be 
separable: 


F(E,L) = L~ 2p f(E), (2) 

where the energy part, f{E), is expressed as jCuddefordll 199ll) 


f(E) = 


2 /3—3/2 


_d_ r B 
d E 


[ (E-*) 


7T 3 / 2 r(m — 1/2 + /3)T(1 — /?) 
, 8 - 3 / 2 +m d m [r 2P p{r)\ 


d<t>" 


2/3—3/2 


/' 


7r 3 / 2 r(m — 1/2 + /3)r(i — p) 

iE-*r 3/ 2 + - dm+ ;r + f )] d$ . 

v ' dT> m + 1 


' max.t) 

Here /? is the velocity anisotropy parameter defined as 


(3) 


g=1 _ {'-'ey - {r e f + {r* 2 ) - (r*) 2 

2((Vr 2 ) - (Vr ) 2 ) 

with v r , ve and v</, being the radial and two tangential components 
of the velocity. The integer, m, is chosen to make the integral con¬ 
verge and depends on the value of /3. In our analysis the parameter 
range of /3 is —0.5 < /3 < 1 and m = 1. /3 > 0 represents radial 
orbits, while tangential orbits have /3 < 0. /3 = 0 corresponds to 
the isotropic velocity distribution. 

In real observations, the tangential velocities of tracers are of¬ 
ten unavailable. We thus test two different cases, in which i) only 
radial velocities are available and ii) both radial and tangential ve¬ 
locities are available. For case (i), the phase-space distribution in 
terms of radius, r, and radial velocity, tv, is given by the integral 
over tangential velocity, vt = \/Vg + 1/2, as 


2 To define the binding energy, we adopt the convention that <b(r) > 0. 


P(r,Vr\C) = J L- 2p f(E)2irvtdvt, (5) 

where C denotes a set of model parameters. With the Laplace trans¬ 
form, this can be written as 


P(r, v r \C) 


I f Er d<fi d r 2p p(r) 

V2Tvr 2 P J #(r- maXit ) \/Er-<& d<E> 


( 6 ) 


where E r = 4>(r) — v 2 /2. All factors of m cancel in the Laplace 
transform and hence Eqn. [6] does not depend on m. For case (ii), 
the distribution function is simply Eqn.(2] i.e. 


P(r,Vr,v t \C) = L~ 2p f(E), (7) 

where E = <I>(r) — v 2 /2 — v 2 /2 and L = rvt- 


3.2 NFW potential and double power-law density profiles of 
the tracer population 

IWilkinson & Evan j j 19991) and lSakamoto et al] d2003l) adopted the 
so-called truncated flat rotation curve model for the underlying dark 
matter pot ential. In our analysis, we wi ll extend Eqn.[2]to the NFW 
potential ( iNavarro'et ~a01996bl H997bh 


i „ 2 ( ln(l + r/r s ) , 

$(r) = —4nGp s r s ( - - - - + 


1 


r/r s 


1 -f r max ,h/r £ 


, ( 8 ) 


when r < r max ,h, and 


<D(r) = —47r Gp s r 


when r > r max ,h- 


2 ( ln(l + r max ,h/r s ) 


x,h/r s 




r/r s 


(r/r s ){l + r max ,h /r s 

(9) 


There are two parameters in Eqn. [8] and Eqn. [9] the scale- 
length, r s , and the scaledensity, p s , defined at r = r s . r max ,h is 
the halo boundary. If the halo is infinite, the second term in Eqn. [8] 
vanishes. In most of our analysis, we will assume the NFW halo 
is infinite. We test different choices of halo boundary in the Ap¬ 
pendix [B] 

To derive analytical expressions for Eqn. [6] and Eqn. [7] we 
need an analytical form for the tracer density profile, p(r). Fig. [2] 
shows the radial density profile of stellar mass (red points) in each 
of the five Aquarius halos. Error bars are obtained from 100 reali¬ 
sations of bootstrap resampling. In most of the cases, these profiles 
can be described well by a double power law (black dashed lines 
are double power-law fits that minimise x 2 )- Significant deviations 
from a double power law are most obvious in the outskirts of the 
halos. For example, halo E has a prominent bump at r ~ 100 kpc 
due to a tidal stream. 

There are indications that the real MW has a two-component 
profile, with density falling off more rapidly beyond ~ 25 kpc, 
whereas M31 has a smooth profile out to 100 kpc with no obvi- 
ous break (e.g. lWatk ins et al.1l2009l:lD eason et all201ll : [Sesar et ahl 
1201 lh . Recently, Deason et alJ 1 2014 ) report evidence for a very 
steep outer halo profile of the MW. If we believe that MW halo stars 
originate from the accretion of dwarf satellites, whether the profile 
is broken or unbroken depends on the d etails of accretion history 
dPeason et al.ll2013l : iLowing et alj|2015h . There is an as yet unre¬ 
solved debate over whether the stellar halo of the MW has an addi¬ 
tional contribution from stars formed in situ, in which case a break 
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in the profile may reflect the transition from in situ-dominated re¬ 
gions to accretion dominated regions. 

As our mock halo stars (which are all accreted) and observed 
MW halo stars can be approximated by a double power-law pro¬ 
file, we adopt the following functional form to model tracer density 
profiles: 


p{r) oc 


— ) + ( — 

roj \r o 


( 10 ) 


This equation has three parameters: the inner slope, a, the outer 
slope, 7 , and the transition radius, ro. 

Previous studies have adopted a single power law to de- 
scribe the densi t y pro fi le of MW halo stars beyond r ~ 20 kpc 
(e.g. IXue et alj [20081: iGnedin et alj l20ld: |5 eason et alj 1201 j. 
Wilki nson & Evan? 1999b . Our double power-law form naturally 
inclu des this p ossibility as a special case. We also note that 
I Sakamoto et al.l (120031) considered the case of “shadow” tracers 
with a radial distribution that shares the same functional form with 
the underlying dark matter. We emphasise that our mock halo stars 
are not “shadow” tracers; their radial distribution is significantly 
different from that of the dark matter. 

Assuming these analytical expressions for <f>(Y) and p(r), 
Eqn.[ 6 ]and Eqn.[7]can be written more explicitly as 


P{r,Vr\p s ,r s ,/3,a,-y,r 0 ) = - 


2/3 ex. 7 n i?max,t ^2^ — 1 


L 


VZnr^Vs ifl inner a/ £ ( ? ') - </>( R ') 

(2/3 - a)(£)“r^ + (2/9 - 7 )&r7 a 


Kf )“»T 7 + (^)-rr“l 2 


-d R' 


( 11 ) 


and 

P(r, v r , v t \ps, r 3 , /3, a, 7 , r 0 ) = 

a - lr 2 P 


2 3/2-/3 7r 3/2 u 3 r ( y0 + l/2)r(l - /3) ' 

/ dR\e{r) - cf>(R')) 0 - 1/2 x 
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( 1 +R’ 

- ln(l + 77')) 
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1 1 
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[(i+fl'H i+fl'J 
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R' 
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_ R ,20+1 _ 

t- — ln(l + R')\ [(f)“rr+(f)%.-“] 
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a 27 } j R‘ 
ro r 0 


J \ CK + 7 —1 


7 2 a ^ j R 1 
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r 0 , 

/ \ CK +7 —1 


+ 


(2/3- 7 )r s -“- 7 

V r 0 r 0 

(2/3 -a)r^-( — ) 2a ' - (20 - 7 K“ 2a - (—' 
ro V ro J ro \ ro J 

( 12 ) 

Here, analogously to Iwilkinson & Evani J 1999h . we have intro¬ 
duced a characteristic velocity, v s = r s y/dnGps. The binding en¬ 
ergy, e, angular momentum, l, potential, cj), and radius, R, have all 


been scaled by v s and r 3 and are thus dimensionless, as follows, 




$ 



(13) 


As mentioned above, i? max ,t is the boundary of the tracer distribu¬ 
tion and, for most of our analysis, we take f? max , t = 00 . Note that 
Eqn.ll 1 land Eqn.[l2]are deduced by assuming the tracer boundary, 
3?max,t, is smaller or equal to the halo boundary, -R max ,h- In both 
Eqn.[TT]and Eqn.[l2]there are six model parameters. 

The phase-space probability of a tracer at radius, r, whose ra¬ 
dial and tangential velocities are v r and Vt, can be derived from 
Eqn.[TT]or Eqn[l 2 ] The lower limit of the integral is determined by 
solving 


</>(i?inner) = t, (14) 

where e equals 4>{R) — v 2 /{2v 2 ) when only the radial velocity is 
available, and e equals <j)(R) — v^/{2v 3 ) — v|/( 2 w 2 ) when tangen¬ 
tial velocity is also available. The fact that the integral goes from 
Rinner to Rmax,t indicates that the phase-space distribution at ra¬ 
dius r has a contribution from tracers currently residing at larger 
radii, whose radial excursion includes r. 


3.3 Likelihood and window function 

The probability of each observed tracer object, labeled i, with ra¬ 
dius, ri, radial velocity, v„ , and tangential velocity, vu, is 

Pi{ri, Vri, vu\p s ,r s , (3, a, 7 , ro). (15) 


Dynamical tracers, such as MW globular clusters, BHB stars 
and satellites, are subject to selection effects. For example, sample 
completeness is often a function of apparent magnitude (hence dis¬ 
tance). If we assume that all selection effects can be described by a 
window function, then the probability of finding each tracer object, 
i, within the data window is given by the normalised phase-space 
density 


/•; = 


p 


P d 3 rd 3 n 


(16) 


The integral in the denominator runs over the phase-space window. 
The likelihood function then has the following form: 


L = HFi. (17) 

i 

It can easily be shown that this likelihood function is equiva¬ 
lent to the extended likelihood function marginalis ed ov er the am¬ 
plitude parameter of the phase-space density (e.g. lBarlowiri990b , 
which we are not interested in. For our mock MW halo star cata¬ 
logue, we deliberately exclude stars in the innermost region of the 
halo. These stars have extremely high phase-space density and so 
make a dominant contribution to the total likelihood, strongly bias¬ 
ing the fit. We find that excluding all stars at r < 7 kpc removes this 
bia0 The window function in our analysis is then simply P = 0 at 
r < 7 kpc. In real observations, the window function can be much 
more complicated. 

We seek parameters that maximise the value of the likelihood 
function defined in Eqn. |T 6 ] and Eqn. |T7] In order to search the 
high-dimensional parameter space efficiently, we use the software 


3 A detailed discussion of the radial dependence of results from our model 
is given in Sec. [6] 
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IMINUIT , whic h is a python interface of the MINUIT function 
minimiser Jjames & Roosfll975l) . 

There are six parameters in Eqn.|TT]or Eqn.[l2] To make best 
use of the likelihood method, we treat all six parameters as free. In 
previous work using this approach the three parameters of the spa¬ 
tial part of the tracer distribution are often fixed to their observed 
values. We have carried out tests and found that, as expected, three 
parameter models give results consistent with those using six pa¬ 
rameters only when the choice of tracer density profile is close to 
the true distribution. We recommend that all six parameters should 
be left free if the observed sample size is large enough to avoid 
introducing unnecessary bias. 

Another source of potential bias in the halo mass estimates 
of previous studies arises from the use of universal mean mass- 
concentration relations for dark matter halos. In CDM simulations, 
the relation between halo mass and concentration has very large 
scatter (e.g. lNeto et al.ll2007h . Taking h alo A as an examp le, if we 
use the mass concentration relation from Duffy et al.1 (120081) . the es¬ 
timated concentration would be around 5.7, which is almost three 
times smaller than the true value (see Table [T}. This would result 
in an overestimate of halo mass by almost an order of magnitude, 
and the corresponding scalelength, r s , would be three times larger. 
The huge scatter in the mass-concentration relation can cause catas¬ 
trophic problems unless we are fortunate enough that the host halo 
of the MW does in fact lie on the mean mass-concentration relation. 

3.4 Weighting tagged particles 

As described in Sec. 12.31 our mock catalogues are created by as¬ 
signing stars from each single age stellar population to the 1% most 
bound dark matter particles in their host halo at the time of their 
formation. The total stellar mass of each population will obviously 
vary from one population to the next (according to our galaxy for¬ 
mation model), as will the number of dark matter particles actually 
tagged (according to the number of particles in the corresponding 
formation halo). The result is that stellar masses associated with 
individual dark matter particles range over several orders of magni¬ 
tude. Particles tagged with larger stellar masses correspond to more 
stars, and thus in principle should carry more weight in the likeli¬ 
hood fit. 

To reflect this we could simply reweight each particle accord¬ 
ing to its associated stellar mass, M*,i. However, individual stars 
are not resolved: the phase-space coordinates of the underlying 
dark matter particles comprise the maximum amount of dynamical 
information available from the tagging technique. Therefore, we 
give each particle a weight This conserves 

the total particle number, 7Vt ags , but re-distributes this among par¬ 
ticles in proportion to the fraction of the total stellar mass they rep¬ 
resent. In this way we maintain a meaningful error estimated from 
the likelihood function. 

We also randomly divide stars into subsamples and apply our 
maximum likelihood analysis to each of these to estimate the ef¬ 
fects of Poisson noise. To do so, we assign each weighted particle 
a new integer weight drawn from a Poisson distribution with mean 
equal to the weight given by the expression above. We repeat this 
procedure 10 times, so that we have 10 different subsamples. The 
expectation values of the total weight for all tagged particles in 
these subsamples are the same, so this approach can be regarded as 
analogous to bootstrap resampling. We find this procedure yields 
consistent error estimates with that obtained from the Hessian ma¬ 
trix of the likelihood surface. From now on, we will only quote 
errors from the Hessian matrix. 



Figure 3. The ratio between input and best-fit halo masses (x-axis) versus 
the ratio between input and best-fit halo concentrations (y-axis). Both ra¬ 
dial and tangential velocities are used. The red cross is the mean ratio over 
750 different realisations, which is very close to 1 on both axes (horizontal 
and vertical dashed lines). Black solid contours mark the region in param¬ 
eter plane enclosing 68.3% (1 <t) and and 95.5% (2 a) of best fit parameters 
among the 750 realisations. 


We restrict our analysis to the 10% oldest tagged particles in 
the main halo. This is to reflect the fact that, in real observations, 
old halo stars such as blue horizontal branch (BHB) and RR Lyrae 
stars are most often used as dynamical tracers, because they are 
approximately standard candles. We also exclude stars bound to 
surviving subhalos. 


3.5 Testing the method 

Before fitting the model distribution function to our realistic mock 
stellar halo catalogues, we test the method with ideal samples of 
particles that obey Eqn. |T2] We applied our maximum likelihood 
method to 750 sets of 1000 phase-space coordinates (r, v r and vt) 
each drawn randomly from the same distribution function of the 
form given by Eqn. G3- ng. a shows a comparison between the 
input halo parameters and the recovered best-fit halo parameters. 
The x axis is the ratio between the best-fit and true-input halo mass, 
and the y axis the ratio between best-fit and true concentration. 
The red cross indicates the mean ratios averaged over all the 750 
realisations, which is very close to unity (horizontal and vertical 
dashed lines). 

The best-fit halo mass and concentration varies among reali¬ 
sations as a result of statistical fluctuations, as shown in Fig. [3] We 
note that the shape of these contours indicate a correlation between 
the recovered halo mass and concentration parameters. The correla¬ 
tion coefficient (i.e., normalised covariance) is -0.89, which implies 
a strong negative correlation in the model parameter. We will dis¬ 
cuss this correlation further in Sec. [5] The above exercise reassures 
us that our method works with ideal tracers, so we can move on to 
apply it to the more realistic mock halo stars in our simulations. 
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Figure 4. The best-fit values of M 200 , C 200 and P for the five halos (black 
dots with errors). Tangential velocities are not used. Error bars are la uncer¬ 
tainties obtained from the Hessian matrix and are almost invisible. The la 
errors are comparable to the scatter among the 10 subsamples constructed 
in Sec. 1531 For direct comparison we show the true values of M 200 , C 200 
and 0 as red dots. 


4 RESULTS 


In this section we investigate how well the true halo mass can be 
recovered by fitting Eqn.QT]to mock halo stars in cases where: (a) 
only radial velocities are available (Sec. 14.1b and (b) both radial 
and tangential velocities are available (Sec. 14.2b . In both cases we 
model the underlying potential with infinite halo boundaries. We 
refer to parameters estimated with the maximum likelihood tech¬ 
nique as best-fit (or measured) parameters, to be compared with the 
real (or true) parameters taken directly from the simulation. 

The total number of tagged particles we used in the five 
halos is shown in Table [2] These are of the order of 10 5 , 
one or t wo orders of magnit ude larger than t he tra cer samples 
used by iDeason et alj d2012h or iKafle et ahi 120121 ). This per¬ 
mits a robust test of the method free from the effects of sam¬ 
pling fluctuations. Future samples of observed tracers will grow 
with ongoing and upcoming surveys such Gaia JPrus tij 1 20121: 
I Gilmore et ah 2012[) and other deep spec troscopic surveys like 
MS-DESI jEisenstein & DESI Collaboration! I2OI5I) and 4MOST 
Ide Jong et a! ■ 2012h. 


Table 2. The total number of tagged particles in each of our five simulated 
halos. 


A 

B 

C 

D 

E 

number 181995 

225030 

184197 

365280 

120806 



r(kpc) 


Figure 5. The velocity anisotropy parameter, 0 , as a function of radius for 
stars (blue solid curve) and a randomly selected subsample of dark matter 
particles (black solid curve) in halo A. Error bars are estimates obtained 
from bootstrap resampling. Blue and black dashed curves are the mean 
value of 0 over the whole radial range for stars and dark matter particles. 


4.1 Radial velocity only 

Fig.0shows, as black points, the best-fit M 200 , C 200 and 0 for our 
five halos in the case where only radial velocities are known. These 
best-fit parameters are given in Table Q] along with the true halo or 
tracer properties (shaded in grey), which are plotted as red points 
in Fig.Q] 

TableQ]lists the true and best-fit values of the host halo mass 
(A/200), halo concentration (C200), scalelength (r s ), scaledensity 
(p s ) and virial radius (/?20o)- Note only two of these parameters 
are independent. The best-fit A/200 values are overestimates of the 
true values for halos B and D by 140% and 7%, and underestimates 
for halos A, C and E by 30%, 10% and 35% respectively. Since the 
number of stars is very large (Tabled, the statistical errors are all 
very small. The systematic biases in the estimates of halo mass and 
concentration are thus very significant. We find the scatter in pa¬ 
rameters among the 10 bootstrap subsamples discussed in Sec. [331 
is comparable to the statistical error in the fit. 

The measured spatial parameters (a, 7 and r 0) agree well with 
the true values obtained from a double power-law fit to the stel¬ 
lar mass density, shown as black dashed lines in Fig. [2] The pro¬ 
files corresponding to the best-fit parameters are plotted as dashed 
green lines in the figure. The agreement is especially good on scales 
smaller than or comparable to the transition radius r 0 . In the out¬ 
skirts, differences become noticeable for halos B and C. This is 
due to the fact that there are fewer stars in these regions and the 
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Table 1. best-fit and true model parameters for each of our five halos. The highlighted rows list the true values of model parameters and the subsequent two 
rows the corresponding best-fit values when using only radial velocities, v r , and using both radial and tangential velocities, v r + vt- 



A 

B 

C 

D 

E 

true M 2 oo[10 12 M q ] 

1.842 

0.819 

1.774 

1.774 

1.185 

V'P 

1.302 ±0.02 

1.972 ± 0.056 

1.616 ±0.029 

1.911 ±0.032 

0.744 ± 0.017 

V r + V t 

1.150 ±0.003 

0.867 ± 0.013 

1.411 ±0.013 

1.410 ±0.005 

0.995 ± 0.014 

true C200 

16.098 

8.161 

12.337 

8.732 

8.667 

V r 

8.616 ± 0.276 

3.080 ± 0.053 

7.682 ±0.140 

6.721 ±0.118 

8.758 ±0.314 

V r + V t 

15.269 ± 0.097 

8.186 ±0.107 

15.878 ±0.199 

10.317 ±0.039 

10.297 ±0.144 

true r s [kpc] 

15.274 

23.000 

19.685 

27.808 

24.493 

v r 

25.422 ± 1.006 

81.660 ± 1.328 

30.642 ± 0.494 

37.035 ± 0.566 

20.752 ± 0.636 

Vr + Vt. 

13.763 ± 0.088 

23.360 ± 0.328 

14.169 ±0.183 

21.805 ±0.086 

19.446 ± 0.287 

true log 10 p s [M 0 /kpc 3 ] 

7.193 

6.591 

7.025 

6.646 

6.642 

Vr 

6.664 ± 0.034 

5.646 ± 0.016 

6.544 ± 0.019 

6.407 ±0.018 

6.681 ±0.038 

V r + V t 

7.278 ± 0.007 

6.610 ± 0.014 

7.321 ±0.014 

6.857 ± 0.004 

6.852 ± 0.015 

true R200 [kpc] 

245.88 

187.70 

242.82 

242.85 

212.28 

V r 

219.025 ± 11.155 

251.525 ±5.962 

235.382 ± 5.786 

248.901 ± 5.786 

181.754 ± 8.579 

Vr ± Vt 

210.144 ± 1.797 

191.234 ±3.383 

224.984 ± 3.785 

224.962 ±1.133 

200.241 ± 3.770 

true /3 

0.660 

0.570 

0.587 

0.752 

0.464 

V r 

0.994 ± 0.001 

1.000 ± 0.001 

1.000 ±0.001 

0.830 ± 0.008 

0.713 ±0.08 

V r ± V t . 

0.458 ± 0.002 

0.397 ± 0.002 

0.407 ± 0.002 

0.553 ± 0.001 

0.254 ± 0.003 

true a 

2.926 

2.912 

3.055 

2.007 

2.223 

v r 

2.965 ± 0.490 

2.911 ±0.007 

3.008 ±0.011 

2.112 ±0.009 

2.454 ± 0.023 

Vr + Vt. 

2.774 ± 0.010 

2.770 ± 0.008 

2.962 ± 0.014 

2.012 ± 0.005 

2.413 ±0.017 

true 7 

6.468 

7.485 

6.383 

6.048 

5.256 

Vr 

6.650 ± 0.037 

8.362 ± 0.038 

5.884 ± 0.034 

6.031 ±0.017 

5.297 ± 0.023 

Vr ± Vt 

6.110 ±0.025 

8.140 ±0.034 

5.623 ± 0.030 

5.820 ± 0.013 

5.305 ± 0.020 

true ro [kpc] 

51.892 

38.506 

60.040 

40.121 

24.215 

v r 

53.376 ± 0.260 

38.779 ±0.138 

57.643 ± 0.847 

42.183 ±0.204 

26.645 ± 0.204 

Vr ± Vt 

42.590 ± 0.345 

35.736 ±0.111 

51.375 ±0.935 

38.165 ±0.100 

26.645 ± 0.294 


profiles have a significant contribution from coherent streams. As a 
result, the direct fitting of radial profiles returns parameters that are 
slightly different from those obtained from the likelihood technique 
because the latter also involves fitting to the velocity distribution. In 
addition, the direct fit is dependent on our choice of radial binning. 

The velocity anisotropy, /?, is poorly estimated. The model as¬ 
sumes f3 to have a single value at all radii. However, the true ve¬ 
locity anisotropy in the simulation does depend on radius: the blue 
solid curve with errors in Fig. [5] is the velocity anisotropy profile 
of stars in halo A as a function of distance from the halo centre. 
We also show the mean value of f3 (0.66 in Table Q} as the blue 
dashed line. The poor measurement of f3 is not simply due to ra¬ 
dial averaging, because we can see that the estimate of /3 for halo 
A, 0.994, is significantly greater than the real value over the whole 
radial range probed. The black curve shows the radial profile of /3 
for dark matter particles. There is an obvious offset between the 
velocity anisotropy profiles of stars (tagged particles ) and all dark 
matter, which we will discuss in detail in Appendix lAl 


4.2 Radial plus tangential velocity 

Best-fit parameters when tangential velocities are also included are 
shown as black points in Fig.[ 6 ]and in Table[I] Compared with the 
results when only radial velocities are used, we see a reduction in 
the overall bias of the best-fit parameters with respect to their true 
values. However, there are still significant discrepancies between 
best-fit and true parameters, compared with the small errors. M 200 
is underestimated for halos A, C, D and E by 40%, 20%, 20% and 
15% respectively. For halo B there is a 5% overestimate. Although 
the measured host halo masses seem to be worse for halos A, C and 


D, compared to the case where only radial velocities were used, the 
agreement between measured and true halo concentrations in the 
same halos is significantly improved. The best-fit spatial parame¬ 
ters, on the other hand, converge to stable values that agree well 
with true values. 

Compared with Fig. [4] the measurements in Fig. [ 6 ] are much 
more clearly correlated with the true halo properties. In particular, 
the shape of the best-fit and true (3 curves are in good agreement, 
although there is approximately a constant offset between them. 
Tangential velocities are therefore essential for measuring tracer 
velocity anisotropy, but even with this information there can still 
be a systematic bias in the absolute value of /3 recovered by distri¬ 
bution function modelling. We return to this point in Appendix lAl 


4.3 Overall model performance 

We have shown that the degree of bias between true and best-fit val¬ 
ues resulting from our fitting procedure differs from halo to halo. In 
the current subsection we aim to show how well the model works in 
recovering the overall phase-space distributions of our mock halo 
stars. Fig. [7] shows phase-space contour plots for mock halo stars 
(green solid curves) and compares them with the predictions of our 
model (red dashed curves). We choose to make the contour plot in 
terms of two observable quantities: kinetic energy, K, and angular 
momentum, L. Each row shows a different halo. In each column, 
the choice of model parameters is different. In the leftmost col¬ 
umn, we use the true values of all parameters. In the case of f3, we 
take its “true” value to be the velocity anisotropy averaged over the 
whole radial range. In the central column, we use the true values 
for all parameters except (3, which is set to be the best-fit value 
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Figure 6. The measured M 200 . C200 and 0- This is similar to Fig. [4] but 
based on both radial and tangential velocities. Black solid dots with errors 
show our best-fit model parameters, while red dots show the true values of 
M 200 , C200 and 0. 


in Table [T| obtained using both radial and tangential velocities. In 
the rightmost column, we set all parameters to their best-fit values 
(again using both radial and tangential velocities). Since the green 
solid curves show data from the simulation, they are identical for 
all three columns of a given row. The contour levels correspond to 
the mass density of the 10, 30, 60 and 90% densest cells. 

The distribution functions defined by the true parameters (red 
dashed curves) are a poor description of the mock stars in the left 
hand column, especially for halos A, C and D where we see a sig¬ 
nificant over prediction of low angular momentum particles. Halo 
E is the exceptional case in which we find good agreement. The 
strongest disagreement in the other halos is, interestingly, mainly 
due to the biased measurement of /3. In the central column, where 
we fix the value of /3 to its best-fit value (obtained using both ve¬ 
locity components) we see that the model predictions agree much 
better with the true phase-space distribution, although some dis¬ 
crepancies remain. 

The fact that the model does not properly represent the dis¬ 
tribution of mock stars when we use the true value of /3 is indica¬ 
tive of possible deficiencies in the model functional form. We will 
show in Appendix A that the physical interpretation of (3 in the 
power-law index, —2/3, of our distribution function as the true ve¬ 
locity anisotropy is not appropriate. Moreover, the approximation 
of a constant f3 over the whole radial range is problematic, as we 


true true + fit/3 best fit 



log 10 lf[(km/s) 2 ] log 10 ff[(km/s) 2 ] log 10 A"[(km/s) 2 ] 


Figure 7. A phase-space contour plot (kinetic energy, I\, versus angular 
momentum, L) of mock halo stars (green solid curves) and model predic¬ 
tions (red dashed curves). The left column is based on true halo parameters, 
true 0 and true spatial parameters of stars. The middle column is identical 
to the left column, except that 0 has been fixed to its best-fit value in Ta- 
ble[T|(when both radial and tangential velocities are used). Halo and tracer 
parameters in the right column have all been fixed to be the best-fit param¬ 
eters with both radial and tangential velocities. Contours for the five halos 
are presented in different rows, as indicated by texts in the left column of 
each row. 

know /3 is radially dependent (see Fig. |5). However, this is very 
likely subdominant because the true value of /3 is above the best-fit 
value (0.458) over almost the entire radial range (blue solid curve 
in Fig. [5}. 

For comparison, the right hand column of Fig. [7] shows that 
model predictions based on the best-fit parameters give an equally 
good match to the simulation data. Judging by eye alone, it would 
be hard to tell whether the middle column shows better or worse 
agreement than the right hand column. However, judging according 
to the likelihood ratios, the best-fit halo parameters are indeed a 
much better description of the data than the true halo parameters 
(2S> 3 <t). This is also reflected in the small formal errors of the fit. 


5 SOURCES OF BIAS 

Fig.[7]mdicates that the model is able to recover the general phase- 
space distribution of the mock halo stars, although there are some 
subtle factors which significantly bias our best-fit parameter values 
relative to their true values. There are several possible sources of 
this bias: 

• Correlations among parameters make the model more sensi¬ 
tive to perturbations and, in some cases, a poor fit to one parameter 
will propagate to affect the others. 

• Tracers may violate the assumption of dynamical equilibrium. 

• Both the underlying potential and the spatial distribution of 
tracers may not satisfy the spherical assumption. 
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Table 3. Correlation matrix of model parameters for halo E 


v r + Vt 


P 

M200 

C 200 

a 

7 

ro 

P 

1.0 

0.025 

0.033 

- 0.00006 

- 0.007 

0.001 

M200 

0.025 

1.0 

- 0.887 

- 0.207 

0.342 

0.040 

C200 

0.033 

- 0.887 

1.0 

0.252 

- 0.338 

- 0.111 

a 

- 0.00006 

- 0.207 

0.252 

1.0 

0.382 

0.768 

7 

- 0.007 

0.342 

- 0.338 

0.068 

1.0 

0.724 

ro 

0.001 

- 0.040 

- 0.111 

0.768 

0.724 

1.0 

v r only 


P 

M200 

C 200 

a 

7 

ro 

P 

1.0 

- 0.267 

- 0.787 

0.573 

0.306 

0.507 

A/200 

- 0.267 

1.0 

- 0.321 

- 0.244 

0.187 

- 0.073 

C200 

- 0.787 

- 0.321 

1.0 

- 0.384 

- 0.415 

- 0.473 

a 

0.573 

- 0.244 

- 0.384 

1.0 

0.574 

0.875 

7 

0.306 

0.187 

- 0.415 

0.574 

1.0 

0.799 

ro 

0.507 

- 0.073 

- 0.473 

0.875 

0.799 

1.0 

• The velocity anisotropy (/ 3 ) is not constant with radius as as¬ 
sumed in the model, although this variation is probably subdomi- 


nant compared with the systematic bias in / 3 . 

• The true dark matter distribution may deviate from the NFW 
model. 

• There are ambiguities in how to model the boundaries of ha¬ 
los. 

We have investigated each of these possibilities and found that 
ambiguity in the treatment of halo boundaries are relatively unim¬ 
portant; hence we describe their effects in Appendix B. We have 
investigated the density profiles of the Aquarius halos and found 
that halo A is not well fit by an NFW profile; instead its inner and 
outer density profiles are better described by two different NFW 
profiles of different mass and concentration. This might explain the 
systematic underestimation of A/200 ■ For the other halos, the NFW 
form is a good approximation and thus deviations from it are not 
the dominant source of bias. 

In Fig. [6] we showed that the velocity anisotropy, / 3 , varies 
strongly with radius. The best-fit value of /3 (which is assumed to 
be constant) turns out to show a significant bias. We will discuss 
the origin of the bias for /3 in Appendix A. In the following, we 
will first discuss whether the bias and approximate treatment of /3 
affects the fit of the other parameters. In the following subsections, 
we focus on correlations among model parameters, the spherical 
assumption and the dynamical state of tracers. 

5.1 Correlations among model parameters 

Fig. [ 3 ] demonstrated a strong correlation between A/200 and 0200- 
From a modelling perspective, this is dangerous: there are multiple 
combinations of halo parameters that can give similarly good fits 
to both the tracer density profile and velocity distribution. In this 
subsection we ask what causes this correlation and whether there 
are similar correlations among other parameters. In particular, we 
have seen that the model gives strongly biased estimates of the ve¬ 
locity anisotropy of stars, / 3 . We want to check whether this bias 
propagates to the other parameters. 

Fig.[8]shows the marginalised 1 , 2 and 3 -er error contours for 
all possible combinations of two model parameters and is for halo 


E (tangential velocities are included as constraints). Fig. [ 9 ]is simi¬ 
lar to Fig. [8] but shows the corresponding error contours when only 
radial velocities are used. All parameters have been scaled by their 
true values (Table Q}. The error contours are obtained by scanning 
likelihood values over the full 6-dimensional parameter space. We 
also overplot as black ellipses the 1 -<t error from the covariance ma¬ 
trix recovered for halo E. The agreement between the black ellipses 
and blue error contours is very good in all the panels, indicating the 
error estimated from the Hessian matrix is robust. The correspond¬ 
ing values of the normalised covariance matrix are also provided in 
Table [ 3 ] 

For the other halos, the error contours look qualitatively sim¬ 
ilar, except for halo A in which the correlation between halo mass 
and concentration is weaker. This is because the bias in the recov¬ 
ered halo properties of halo A is mostly due to its deviation from 
an NFW profile. Otherwise, we found the agreement between the 
error ellipse from the covariance matrix and the scanned error con¬ 
tours are worse for the radial velocity only case of halo A, B and 
C. This is mainly because the best fit value of /3 touches the upper 
boundary /3 = 1 and thus the quadratic approximation is no longer 
good enough to estimate the error from Hessian matrix, hut even in 
such cases the errors estimated from Hessian matrix is still accept¬ 
able with at most a factor of two deviation from the more strictly 
obtained error contours. 

From Fig. [8] Fig. [ 9 ] and Table [ 3 ] we can see the correlation 
between A/200 and C200 is very strong when including tangential 
velocities (covariance close to - 1 ). The correlation is not as strong 
if only radial velocities are used. To help understand this correla¬ 
tion, we have explored the velocity distribution of tracers predicted 
by the model using different combinations of A/200 and 0200- We 
verify that, if M200 is increased, the predicted velocity distribution 
of stars extends to larger velocities, with a corresponding reduc¬ 
tion in the probability of smaller velocities. A decrease in C200 can 
roughly compensate for this change in the velocity distribution. 

It is worth emphasising that although the error contours for 
M200 and C200 are highly elongated (corresponding to the correla¬ 
tion between the two), they are still closed, indicating the constrain¬ 
ing power is not insignificant. Because these contours represent the 
statistical error, they can be reduced by increasing the sample size. 
With our current sample of 10 s particles, the statistical errors on 
M200 and C200 are controlled to the 1% level, which is negligible 
compared to the systematic bias in the parameter values. In other 
words, the true M200 and C200 values lie well outside the 3 <r con¬ 
fidence contour, so that statistical fluctuations do not explain the 
deviation between the fit and the true parameters even after consid¬ 
ering the correlation. 

It is interesting to see that the best fit A/200 and C200 in Fig. 6 
tend to be biased in opposite directions, except for halo A. Such bi¬ 
ases are mainly systematic, because the statistical errors are much 
smaller. This indicates a negative correlation between the system¬ 
atic biases, similar to the statistical correlation we have seen above. 
Note that in principle the correlation in the systematic biases could 
happen along any direction, independent of the statistical correla¬ 
tion, and it is unclear why the two act in the same direction here. A 
clean and thorough exploration of this involves segregating various 
model assumptions and adopting a large sample of halos; Hart et al. 
d2015allbh present part of this work, which will also be discussed in 
a forthcoming paper by Wang et al. (in preparation). At this stage, 
we provide some further discussions on this point in our conclu¬ 
sion. 

Similar correlation between the velocity anisotropy parameter, 
/ 3 , and halo properties have been discussed in some previous stud- 
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Figure 8. The marginalised error contours for all possible combinations of every two out of the six model parameters (halo E). Both radial and tangential 
velocities are used. Red, blue and purple contours correspond to 1, 2 and 3-cr errors. Green dots show the best fit parameters scaled by their true values. The 
1-cr error ellipses estimated from the covariance matrix are over-plotted in black. 


ies of the Milky Way and dwarf galaxie s (e.g. IWalker et al.ll200 l j : 
IWolf et all 1 20 id ; iNesti & S aluccil 20131). th ough their models are 
different. In particular. IWalker et all d2009h and IWolf et all d201Gh 
have reported the mass within the half light radius of dwarf galax¬ 
ies is relatively insensitive to the value of j3. Here we have also 
tested whether our model can better constrain the mass within the 


half-mass radiufl of our stellar tracers, and the results are shown in 

Fig.uoi 

The two panels of Fig. [TO] show results based on both radial 
and tangential velocities (top) and radial velocities only (bottom). 
Black dots with errors are the ratio between the best fit mass and 
the true mass within the half mass radius. Encouragingly, we see a 
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4 The half-mass radius is defined to be the radius inside which the enclosed 
stellar mass is half of that of the whole tracer population. 
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very good agreement between the best fit and true mass if tangential 
velocities are used. The levels of biases are about 3.8%, 0.7% and 
2.4% underestimates for halos A, C and D. For halos B and E the 
mass is overestimated by 0.1% and 0.5% respectively. On the other 
hand, if only radial velocities are used the bias is much more signif¬ 
icant. We underestimate the mass by about 27.8%, 33.2%, 31.4%, 
11.1% and 23.8% for the five halos. Compared with Fig.[4] the level 
of bias becomes significantly smaller for halo B, and is slightly im¬ 
proved for halos A, C and D. Our results suggest if tangential ve¬ 
locities are available, the mass within a fixed radius close to the 


half-mass radius of tracers are not sensitive to the parameter corre¬ 
lation and can be constrained much better than the total halo mass. 

In Fig. El we examine the halo mass profiles (cumulative 
mass within a certain radius as a function of the radius) with the 
best-fit parameters (when both radial and tangential velocities are 
included), normalised by profiles with the true parameters. Except 
for halo B. which gives an acceptable result at all radii, the mea¬ 
surements are very close to the true value for r ^ O.2.R200 with a 
less than 5% bias, though the bias is still significant given the small 
statistical errors. The measurements become significantly more bi¬ 
ased at larger radii. The vertical lines mark the locations of half- 
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Figure 10. The ratio between the best fit and true mass enclosed within 
half-mass radius. Errors are obtained through error propagation from the 
covariance matrix of p 3 and r 3 , with the correlated error between p a and 
r s included. 



Figure 11. The best-fit total mass within a fixed radius compared to the 
true mass within that radius. Both radial and tangential velocities are used. 
Errors are obtained through error propagation from the covariance matrix of 
p B and r s , with the correlated error between p B and r a included. The black 
dashed line marks equality between the measured and true mass. 


mass radii of stellar tracers in the five halos, which are close to 
0. 1 -R 200 ■ In practice this means the mass enclosed within 40 to 60 
kpc can be more robustly determined as it suffers less from the 
correlation, and this is usually the radial range for which relatively 
large numbers of tracers can be observed. Thus the mass measure¬ 
ments reported within 40 to 60 kpc in the literature are expected 
to be more robust. Note, however, the result here is obtained from 
stars over the whole radial range. We will further show in Sec.[ 6 ]that 
the mass within about 0 . 2 f ?200 can also be determined robustly if 
stellar tracers are restricted to be within 60 kpc. 

In contrast to the strong correlation between M 200 and C200, 
the correlation between j3 and all the other parameters is very weak 
if tangential velocities are included. This is fortunate, as it suggests 
the systematically biased estimate of j3 will not introduce further 
bias to the other parameters if we have the proper motion informa¬ 
tion. On the other hand and for the radial velocity only case, the 
correlation between (3 and the other parameters are actually quite 
strong, especially in terms of the correlation with halo concentra¬ 
tion. This suggests if one only has radial velocity information, it is 
hard to get a robust estimate of /3 (Fig. [4} and the bias in j3 may 
affect the fitting of the other model parameters. Compared with the 
systematic bias in f3, the radial dependence of /3 is actually sub¬ 
dominant. 

Correlations between halo parameters (M 200 or C200) and 
tracer spatial parameters (a, 7 and r 0 ) are at the level of a few tens 
of percent. For the case when tangential velocities are included, an 
increase in the tracer density outer slope would cause an increase in 
the recovered halo mass and a corresponding decrease in halo con¬ 
centration; conversely, an increase in the inner slope would cause 
a decrease in the halo mass and an increase in the concentration. 
As a result, uncertainties in the fit to the tracer density profile may 
further bias the best fit halo parameters. For example, the best-fit 
(green dashed) curve in the halo C panel of Fig.[2]agrees well with 
the true profile (red points) inside 170 kpc but is somewhat shal¬ 
lower at larger radii. If we fix the three spatial parameters in our 
fit to halo C to those given by a conventional reduced-y 2 best-fit 
to the tracer density (dashed black curve in Fig. [2} the best-fit halo 
mass is boosted by about 10%. If the tracer density profiles deviate 
from the double power-law form, these correlations between halo 
parameters and spatial parameters would introduce further bias to 
the best-fit halo mass. 

Lastly, we note that correlations between the three spatial pa¬ 
rameters are strong as well. This quantifies our earlier finding that, 
in the case of halo A, adding tangential velocities as constraints 
in the fit makes the outer slope of the tracer density profile shal¬ 
lower and the break radius smaller, but results in very little per¬ 
ceptible change in the overall profile shape. Hence, good fits to the 
tracer density profiles may not be unique. An increase in r 0 can be 
roughly compensated by a corresponding increase of both a and 7. 

5.2 Model uncertainties in the spherical assumption 

In our analysis and the majority of existing studies of using dynami¬ 
cal tracers to constrain the MW host halo mass, both the tracer pop¬ 
ulation and the underlying potential are modeled assuming spher¬ 
ical symmetry. However, we know dark matter halos in N-body 
simulations are triaxial (e.g. ljing & Sutol2002l) . and the spatial dis¬ 
tribution of tracers is unlikely to be perfectly spherical either. It is 
thus necessary to investigate how the triaxial nature of the underly¬ 
ing dark matter potential and tracer populations affect our results. 

To do this, we first rotate the five halos to a new Cartesian 
coordinate system defined by their principle axes. In this rotated 
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Figure 12. The host halo masses obtained by fitting to samples of stars 
drawn from different sky directions and using both radial and tangential 
velocities. Halos have been rotated based on their dark matter distribution. 
The £-axis and 2 -axis in the rotated system are chosen to be the major and 
minor principles of the dark matter halo. Measurements are repeated for six 
survey cones centred at =brr, =by and ±z direction, with opening angle d= J. 
Results based on all stars are plotted as the right most point. The red dashed 
lines show the true values of M 200 • 


system, the 2 -axis is aligned with the minor axis of the halo and 
the x-axis with major axis of the halo. We then repeat our analysis 
using six subsamples of tracers drawn from mock ‘survey’ cones 
pointing along each of the three axes from the origin at the centre 
of the halo, in the positive and negative directions. The opening 
angle of each cone is d= j . 

Fig. CO shows the recovered host halo masses for each of the 
six cones. Tangential velocities are included in the fit. For a direct 
comparison, we have also plotted results based on all tracers, as the 
right most point (from TableQ}. There are significant variations in 
the results obtained from surveys along different directions, rang¬ 
ing from only a few percent (halo A) to as much as a factor of two 
(halo D and E). Halos D and E have the most obvious variations. 
We have explicitly checked that the significant overestimate along 
the positive y -axis of halo D is due to the existence of four rel¬ 
atively massive subhalos (M su bhaio/A/ 2 oo,host > 1 %), while the 
large variation in halo E is due to one prominent stream (see the 
yellow dots in the bottom panel of Fig.Q3]or Fig. 6 in lCooper et al.1 
OOlOh ). which ranges from r ~ 80 kpc (~ 0.3f?20o) all the way to 
the virial radius. 

These variations are, however, almost random and uncorre¬ 
lated with the choice of any particular principle axis, and they 
change from halo to halo. Halo A has the smallest variation, with 
all results well below the true host halo mass (red dashed line). Al¬ 
though stronger variations are seen for halo C, all results are again 
well below the true host halo mass. Thus despite the fact that the 
variations between directions can be as large as a factor of two, this 
does not seem to be the dominant cause of the systematic differ¬ 
ences between the best-fit and true halo mass in our analysis. 


Figure 15. Similar to Fig. [TJ but only using stars whose parent satellites 
have been entirely disrupted. 

5.3 Unrelaxed dynamical structures 

The model distribution function used in our analysis assumes that 
the tracer population is in dynamical equilibrium and hence the 
phase-space density of tracers is conserved. Our mock halo stars are 
all accreted from satellite galaxies, with a range of accretion times. 
Some prominent phase-space structures, such as stellar streams, 
may therefore violate the assumption of dynamical equilibrium. In 
this section we ask how the presence of unrelaxed dynamical struc¬ 
tures affects our results. 

We expect the dynamical state of stars in our mock catalogue 
to depend on the infall redshift of their parent satellite, at least ap¬ 
proximately (satellites on different orbits will have different rates 
of stellar stripping). We might expect to obtain an improved mass 
estimate if we use only stars from satellites that fell in earlier, be¬ 
cause they have had more time to relax in the host potential. To test 
this, we rank halo stars according to their infall timcQ We measure 
the host halo mass and concentration with samples defined by dif¬ 
ferent cuts in infall time using both radial and tangential velocities, 
corresponding to roughly the same fraction of stellar mass in each 
halo. The top and middle rows of Fig. |T3] present these parame¬ 
ters as a function of the fraction of stars selected by each cut, for 
the five different halos. A small fraction corresponds to an earlier 
mean infall time, but also (obviously) to a smaller sample size. A 
fraction of 1 means all the mock stars have been included, hence 
the corresponding parameters are those listed in TableQ] 

We see fluctuations in the measured halo properties with in¬ 
fall time, but no obvious trends. Using samples of stars with earlier 
mean infall times does not seem to reduce the bias between best- 
fit and true parameters. This may be because the dynamical state 
of tracers depend on many other factors, such as the orbit of their 

5 Defined as the simulation output redshift at which the parent satellite of 
each star reaches its maximum stellar mass, which is generally within one 
or two outputs of infall as defined by SUBFIND. 
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Figure 13. The host halo masses (M200) and concentrations (C200) measured using a certain fraction of stars that have the earliest infall time and with both 
radial and tangential velocities. The five columns are for five Aquarius halos, as labeled at the top. The red dashed lines show the true values of A/200 and C200 ■ 
The top and middle rows are based on the 10% oldest tagged particles. The bottom row is analogous to the top one, except that stars whose parent satellites 
have not been entirely disrupted yet are excluded. 


parent satellites(El Samples for which the measured halo masses in¬ 
crease produce a corresponding decrease in the measured concen¬ 
trations, again reflecting the strong correlation between A/200 and 
C200- 

To gain more intuition regarding the dynamical state of halo 
tracers, Fig.ll4lshows phase-space scatter diagrams for mock halo 
stars (radius, r, versus radial velocity, v r ). Points are colour coded 
according to the infall time of their parent satellite, with black 
points corresponding to satellites falling in earliest and blue, ma¬ 
genta, red and yellow points to successively later infall times. Stars 
with earlier infall times are clearly more centrally concentrated. For 
points in Fig.[ 6 ]with decreasing fraction of stars that fell in earliest, 
the corresponding particles in Fig. M can be found by excluding 
yellow, red, magenta and blue points by sequence and looking at 
the remaining points. 

Green curves in Fig. |T4] are contours of constant angular mo¬ 
mentum and binding energy. There are six contours in total, cor¬ 
responding to three discrete values of binding energy and two dis¬ 
crete values of angular momentum: dashed lines have a higher an- 

6 We have carried out an analogous exercise in which we rank stars by the 
time at which they are shipped from their parent satellite. We found that 
this stripping time correlates with the infall time of the parent satellite, and 
the conclusions regarding the recovered halo parameters are similar. 


gular momentum than solid lines. Smaller maximum radius indi¬ 
cates higher binding energy. It is thus straightforward to see that 
particles with higher binding energy have smaller velocities and are 
more likely to be found in the inner regions of the halo. Compar¬ 
ing the solid and dashed contours, we see that increasing angular 
momentum at fixed binding energy causes significant differences 
in the inner regions of the halo, while at larger radii the two sets of 
contours almost overlap. 

We can see that points with the same colour trace these con¬ 
tours with some scatter, implying that stars whose parent satel¬ 
lites fall in at a particular epoch share similar orbits. This can be 
seen more clearly in the bottom right panel, which shows a scatter 
plot of binding energy versus angular momentum for stars in halo 
A. Points with the same colour occupy regions covering a narrow 
range of binding energy. The correlation between infall time and 
binding energy of subhalos has been studied bv lRocha et alJd 2012 h . 
Here we have shown that stars from stripped subhalos show a cor¬ 
relation between infall redshift and binding energy as well. 

Although mock stars trace the green contours overall, we can 
still see some prominent structures. For example, there are two yel¬ 
low spurs in the outskirts of halo D and one yellow spur in halo 
E. These correspond to particles that have only just been stripped 
from their parent satellites. These stars are far from equilibrium: 
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Figure 14. Phase-space scatter plot (radial velocities versus radial positions) for the 10% oldest tagged particles in the five Aquarius halos. Data points are 
colour-coded according to their infall time. In sequence, points with yellow, red, magenta, blue and black colours are stars that have earlier and earlier infall 
times. The bottom right panel is a scatter plot of energy and angular momentum for halo A. Only one out of every 20 points are plotted in order to avoid 
saturation. Green curves are contours of constant angular momentum and binding energy (see text for details). 


their exclusion causes the rapid change in M 200 and C 200 in Fig. ED 
between fractions of 1 and ~ 0.7 in halos D and E. 

To confirm that these unrelaxed phase-space structures can af¬ 
fect our results, we have repeated the above exercise excluding all 
stars whose parent satellites have not been entirely disrupted. Cor¬ 
responding results are shown in the bottom row of Fig. ED again 
ranking stars by their infall time. Measured halo masses are clearly 
affected by excluding stars whose parent satellites still survive. For 
halos A and C, we see some small fluctuations in the measured 
halo mass, but the systematic underestimate of the true halo mass 
still remains. The most dramatic changes occur for halos B, D and 
E. First, the point corresponding to a fraction of 1 for halos D show 
a significant increase in the recovered mass towards the true val¬ 
ues, reinforcing our conclusion that unrelaxed structures are caus¬ 
ing significant underestimates of M 200 in these halos. In fact, most 


of the yellow dots in halo D panel of Fig.[l4]are stars that have been 
stripped from satellites that still survive. After excluding these, the 
two highest-fraction points in the halo D panel are based on almost 
the same sample of stars. We also notice that fluctuations around 
the true value for the different fractions are reduced in the bottom 
row (for example, the two lowest fraction points in halo D). 

The effects of excluding halo stars from surviving satellites are 
more ambiguous in halos B and E. The recovered mass of halo B 
decreases slightly, while for halo E the rightmost point, correspond¬ 
ing to all stars, increases, but the two left most points decrease in 
amplitude, causing a stronger deviation away from the true values. 

We further investigate how the measurements pointing in dif¬ 
ferent directions change with more relaxed stars. Fig. Q3] repeats 
Fig. ED using only those stars from satellites that have been en¬ 
tirely disrupted. The measured M 200 of halo A shows some small 
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variations compared with Fig. ED but the variation is too small to be 
significant. The recovered mass of halos B and C improves in some 
directions, whereas in some other directions it worsens. The most 
encouraging improvements are for halos D and E. The two mea¬ 
surements along ±at directions of halo D remain almost unchanged, 
while the measurements in all the other directions are significantly 
improved. For halo E the recovered mass increases significantly in 
all directions. 

Our conclusion is thus for halo D (and perhaps E) the un¬ 
derestimates of their host halo masses when all particles are used 
are mainly due to unrelaxed dynamical structures; for the other ha¬ 
los, the effects of unrelaxed dynamical structures are not obvious. 
Stars with surviving parent satellites in halos A, B and C could 
be more dynamically relaxed and thus excluding stars that are ex¬ 
pected to be unrelaxed does not make a signifi cant change. Further¬ 
more, we note in a recent work. lBonaca et ali (120140 has developed 
a new method of determining halo potential using tidal streams. 
They found individual streams can both under- and overestimate 
the mass, but the whole stream population is essentially unbiased. 
Though their method is different from ours, it is possible that a 
single dynamically hot stream can potentially bias the result, the 
combination of several streams can help to bring an overall unbi¬ 
ased measurement. A more detailed stud y quantifyin g the dynami¬ 
cal state of tracers has been carried out in lHan etail 1 2015al lbti. 



Figure 17. best-fit host halo masses as a function of the outer radius limit of 
the tracer population (r < r cu t). Both radial and tangential velocities are 
used. 


6 MODEL UNCERTAINTIES IN THE RADIAL AVERAGE 
AND IMPLICATIONS FOR REAL SURVEYS 

We have seen in the previous section that our maximum likelihood 
technique recovers different halo mass from sets of tracers with 
different infall redshifts, or more fundamentally, different binding 
energies. The sense and magnitude of these differences show no 
obvious correlations with either quantity, however. Stars falling in 



Figure 16. Best-fit host halo masses using samples of stars in four radial 
bins, (7-20) kpc, (20-50) kpc, (50-100) kpc and >100 kpc. Both radial and 
tangential velocities are used. 


earlier typically have high binding energy and are mostly concen¬ 
trated in the central regions of the halo; since binding energy corre¬ 
lates with radius, we may also expect fluctuations in the recovered 
halo parameters when using samples of stars drawn from a partic¬ 
ular radial range. In this section we investigate this radial depen¬ 
dence. This helps to understand the behaviour of the full model, 
which averages over all radii. Variations with radius are relevant 
to observational applications as well, because in practice tracers 
are often selected from relatively narrow radial ranges, and these 
ranges may be different for different tracers. 

We assign stars to four bins of galactocentric radius: (7- 
20) kpc, (20-50) kpc, (50-100) kpc and >100 kpc. Note in our 
measurements stars inside 7 kpc have been excluded (Sec. 13.3b . The 
model distribution function is then fit to stars in each bin separately. 
However, in each case the three spatial parameters of the tracer dis¬ 
tribution are fixed to their best-fit values obtained from tracers over 
the entire radial range, otherwise we would end up with extremely 
poor extrapolations based on the local density slope. All the other 
parameters, M 200 . C200 and /3, are left as free parameters. The win¬ 
dow function, Eqn.[l 6 ] is modified appropriately for each bin. 

Fig. |T6] shows the measured M 200 as functions of the mean 
radius of each bin, normalised by the halo virial radius (Rioo)- The 
value of M 200 varies significantly with the tracer radius. Thanks to 
the large number of stars (Table |2jl, the errors are all quite small in 
each bin. For halos A, C and E, stars in the outermost (r > 100 kpc) 
and innermost (r < 20 kpc) bins give underestimates, while stars 
at 20< r <100 kpc give significant overestimates. Similarly, for 
halos B and D, stars at r >100 kpc give underestimates, whereas 
stars at smaller radii give overestimates. 

The velocity anisotropy of tracers, /?, is a function of radius, 
whereas the model distribution function assumes a single value of 
/3. To test whether the radial average of fi may affect our estimates 
in the host halo mass, we repeated the analysis of Fig.[l 6 ]but fixed 
the value of /3 in each radial bin to either the best-fit value or the true 
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value obtained from the whole population. These measurements are 
almost identical to the measurements presented in Fig. [16] which 
confirms our result from Sec. ED that the radial averaging of /? 
does not cause further bias in the other parameters. 

One feature in Fig. [T6] is puzzling at first glance: the best-fit 
halo masses obtained from the four radial bins individually are all 
larger than the best-fit halo mass (A /200 = 1.15) obtained using 
stars over the whole radial range. This seems odd, as we might 
expect that the best-fit A /200 would be close to the average of the 
values estimated from the four radial subsamples. The situation is 
not that straightforward, however, because the likelihood surfaces 
from the subsamples are superimposed in two dimensional ( M 200 
and C 200 ) space when the full sample is used. Coupled with the 
strong correlation between the two parameters, the peak of the final 
likelihood surface is located around a region where the correlation 
lines from different subsamples intersect. 

In real observations, there is often a maximum radius of trac¬ 
ers corresponding to the instrumental flux limit. In the present lit¬ 
erature this limit is typically much smaller than the expected halo 
virial radius. Beyond this maximum radius, extrapolations are re¬ 
quired to fit the distributions of both the dark matter and tracers. 
We explore the implications of this directly in Fig. [IT] by adopt¬ 
ing several outer radial cuts (r < r C ut) and reporting the estimated 
halo mass as a function of the cut radius normalised to the virial 
radius (r cu t/i? 20 o)- Unlike Fig.[T 6 ] the three spatial parameters are 
treated as unknown and left free to be constrained by the fit, in or¬ 
der to mimic real observations where the density profiles of tracers 
is taken directly from the available data. 

The overall trends with r cu t in Fig. [IT] are very clear: the re¬ 
covered halo mass is constant at large r cu t, and turns up once r cut 
becomes small (about 40 per cent of i? 20 o)- We checked the best-fit 
values of the tracer spatial parameters in each case, and found they 
do not vary much with the radial cut as long as r cu t < 0.4/?2oo- 
This is because the break radius of tracer density profiles in our 
mock catalogues are smaller than OAR 200 for all the five halos, 
and so the extrapolation in tracer density is not severe. However, 
once r cu t reduces below OAR 200 , the outer slope becomes essen¬ 
tially unconstrained. We believe the turn-up behaviour is due to the 
changing dynamical state of tracers and the extrapolations required 
to know the underlying potential where there are no tracers. 

Previous constraints on the MW halo mass have been de¬ 
rived from tracers roug hly cove ring the range 0.1 to 0.4 R 200 
(R 200 ~250 kpc: iDeason et alJl2012h . Our results suggest that 
the halo mass, A/ 200 , derived from the fitting distribution function 
of these tracers may be significantly biased even with respect to 
‘asymptotic’ results from the same method using all stars in the 
halo. Furthermore, instead of being a sharp cut, the radial selection 
functions of real surveys are often complicated, with non-trivial in¬ 
completeness as functions of radius and angular position. These se¬ 
lection effects may cause additional bias in the measured host halo 
mass. 

We have shown in Fig. II1 1 that the total mass within the half¬ 
mass radius of the stellar tracer population can be constrained more 
precisely than the total mass of the halo. We now test if this conclu¬ 
sion is robust to changes in the radial range spanned by the stellar 
tracers. We repeated the measurements shown in Fig.[TT]using only 
the stars within 60 kpc. The results are displayed in Fig. [T8] in the 
case when both radial and tangential velocities are included in the 
analysis. The total inferred mass within a fixed radius is strongly 
biased if this radius is close to the virial radius, R 200 , but, encour¬ 
agingly, as the radius is decreased, the measured enclosed mass 
becomes increasingly close to the true value. Our conclusion that 



Figure 18. The best-fit total mass within a fixed radius compared to the 
true mass within that radius, which is obtained from tracers inside 60 kpc 
and with both radial and tangential velocities. The black dashed line marks 
equality between the measured and true mass. Vertical lines mark the loca¬ 
tion of half-mass radii. 

the mass enclosed within the half-mass radius can be constrained 
reliably still holds even when stellar tracers inside only 60 kpc are 
considered. 


7 CONCLUSIONS 

Several authors have measured parameters of the host halo of 
our MW, in particular its total mass, by fitting specific forms 
of the distribution function to the observed distances and veloc¬ 
ities of dynamical tracers such as old BHB and RR Lyrae stars , 
globular clusters and satellite galaxies dWil kinson & Evans| [l99^ : 
ISakamoto et al.ll2003l : IDeason et al Jl 20 1 3 ) . These models assume 
that the tracers are in dynamical equilibrium within the host poten¬ 
tial. With the help of Jeans theorem, the distribution function of the 
tracers is further assumed to depend only on two integrals of mo¬ 
tion, the binding energy, E, and the angular momentum, L. In the 
case of a separable function of E and L, the distribution function 
can be obtained through Eddington inversion of the tracer density 
profile. 

In this paper we have extended earlier analytical forms of the 
MW halo distribution function to the case of the NFW potential, 
which is of most relevance to CDM-based models. We generalised 
the radial distribution of tracers (halo stars) to a double power law, 
which is suggested by recent observational results and simulations. 
We used a maximum likelihood approach to fit this model distri¬ 
bution function to a realistic mock stellar halo catalogue of dis¬ 
tances and radial velocities, constructed from the high resolution 
Aquari us TV -body simulations using the particle tagging technique 
of lCooper et alJ < 120101) . Our aim was to test the model performance 
and assumptions. We considered cases with and without additional 
tangential velocity data. Our conclusions are as follows: 
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• The best-fit host halo virial masses and concentrations are bi¬ 
ased from the true values, with the level of bias varying from halo 
to halo. 

• Adding tangential velocity data substantially reduces this bias, 
but does not eliminate it. For example, for halo B the agreement 
between measured and true halo mass is very good (a 5% overesti¬ 
mate) if tangential velocities are used, but for halo A, a 40% under¬ 
estimate persists even with this additional constraint. The inclusion 
of tangential velocities therefore is crucial for accurate measure¬ 
ments of both host halo and tracer properties, especially for the 
velocity anisotropies of the tracers. 

• A strong negative correlation between the host halo mass and 
the halo concentration is found in our analysis. 

• The model gives a strongly biased measurement of the veloc¬ 
ity anisotropies of stars. 

• If tangential velocities are available, the correlation between /3 
and all the other parameters are very weak. If only radial velocities 
are used, p is strongly correlated with other parameters and the 
bias in /3 will be propagated to these parameters. This is because 
when tangential velocities are not available, we have to rely on the 
model functional form to infer the unknown tangential component 
and hence /?. 

• Various sources contribute to the biased estimates of halo 
properties. Violation of the spherical assumption is relatively sub¬ 
dominant for the five Aquarius haloes. Violation of the dynamical 
equilibrium assumption, caused for example by streams, could af¬ 
fect the fits significantly, although we do not observe a systematic 
sign for the bias in M200 (that is, unrelaxed substructures cause 
underestimates in some cases and overestimates in others). 

• When including tangential velocities, the systematic bias 
tends to happen along the correlation direction of M200 and C200 
except for halo A. 

• In contrast to the significantly biased measurements of M200 
or C200, the model gives good constraints on the total mass within 
the half-mass radius of stellar tracers when including tangential ve¬ 
locities. 

The strong correlation between M 200 and C200 arises because 
changes in the corresponding tracer velocity distribution due to the 
increase of one of these parameters can be roughly compensated by 
the other. The correlation between M200 and C200 is not as strong 
for the radial velocity only case, which is probably overwhelmed by 
the strong correlation between /3 and the other parameters, reflect¬ 
ing the fact that the dominant source of bias is the model dependent 
fit of the tangential component. If the model fails to properly reflect 
the true phase-space distribution of tracer objects, the best fit P and 
other parameters will be strongly biased. 

There are different combinations of halo parameters which 
give similar likelihood values along the correlation direction. Thus 
the model is vulnerable to perturbations (for example from dynami¬ 
cally hot structures). This can be seen from Fig. [6] the best fit M 200 
and C200 are offset in opposite directions with respect to their true 
parameter values. This is not the case for halo A, because the domi¬ 
nant source of bias for halo A is the deviation from the NFW model, 
and the error contour of M 200 and C200 is not as elongated as the 
other halos. More detai led dis cussions about this halo will be given 
in a future study dHan et ali2015bl) . 

It is, however, confusing to see that although the systematic 
bias tends to happen along the correlation direction of M200 and 
C 20 O) it is much larger than the statistical errors. For example, we 
can see clearly in Fig [8] that the best fit M 200 is about 15<r away 
from its true value. This is probably because the statistical error 


in our analysis does not account for the correlations introduced by 
phase-space structures or clumps. Our mock stellar halo catalogue 
contains a very large number of stars. However, these individual 
stars are not completely independent of each other. Structures such 
as coherent streams are highly correlated in phase space and it 
is possible the true number of independent components is much 
smaller than the total number of stars. To quantify the true number 
of degrees of freedom by considering correlated phase structures 
is beyond the scope of our current study. A more detailed discus¬ 
sion is given by Han et alj l i2015al lbh. in which we introduce a new 
method based on the steady state assumption, independent of any 
other assumptions about the model functional form. 

Encouragingly, we found the mass within the half-mass radius 
of the tracer population to be relatively insensitive to the parameter 
correlations and can be constrained more robustly once tangential 
velocities are used. This is true even when only stars within about 
60 kpc are available. Similar correlations between model parame¬ 
ters and the robustness of the best constrained mass within a fixed 
radius have been reported and discussed in previous work (e.g. 
iDee & Widrowl20l4lKafle et al.l2014l : IWolf et all20l(il) . although 
these models are quite different from ours. For our model, the cor¬ 
relation could be closely related to the fact that there are relatively 
few stars outside O. 2 R 200 , beyond which the stellar radial profiles 
drop very quickly. Given the large number of dynamical tracers in¬ 
side O. 2 R 200 , it is not surprising to find that the mass within this 
radius can be better constrained. In contrast, the total halo mass, 
M 200 , is dominated by mass in the outskirts of the halo and more 
tracers at large radii are required to have better constraints. 

Further information needs to be incorporated into the model 
to weaken the correlation between mass and concentration. For 
example, including more tracers at large radii (perhaps from tidal 
streams) may help to weaken the correlation and improve the mea¬ 
surements of mass in the outer halo. Satellite galaxies and glob¬ 
ular clusters in the outer parts of the halo with proper motion 
measurements could be useful. Having two populations of trac- 
ers at two different radial ranges could also be very helpful (e.g. 
[Walker & Penarrubia|[201 ll ). This would enable us to constrain the 
mass at two different half-mass radii and hence the entire mass pro¬ 
file can be fixed. More detailed investigations regarding the nature 
of c orrelation between M 200 and C200 have been carried out by 
Ha n et a l ■ l 2015allbl) . 

The correlations between velocity anisotropy, /3, and all other 
parameters are very weak when including tangential velocities. 
This is fortunate, because we know that the model can give sys¬ 
tematically biased estimates of /3 for stars; this particular bias is 
not propagated to the other parameters when including tangential 
velocities. However, if only radial velocities are available, the con¬ 
dition becomes very different. The correlation between P and all 
the other parameters is strong. Combined with the biased measure¬ 
ments of P in Fig. [4] this suggests that if proper motions are not 
available, it will be difficult to obtain robust constraints on p and 
the bias may affect the fitting of the other parameters. Only by in¬ 
cluding tangential velocities can the correlation between P and the 
other parameters be broken and P be better constrained. 

In addition to the correlation between halo mass and concen¬ 
tration, relatively weak but still significant correlations exist be¬ 
tween these halo parameters and the three parameters describing 
the spatial variation of tracer density. When including tangential 
velocities, we found that a steeper inner slope gives a lower es¬ 
timate of halo mass, while a steeper outer slope gives an higher 
estimate. If the true tracer density profile deviates from the double 
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power-law form, the resulting bias will be propagated to the best-fit 
values of the halo parameters. 

The model distribution function requires tracers to be in dy¬ 
namical equilibrium, with time-independent phase-space density. 
In reality, stars stripped from satellite galaxies can have highly cor¬ 
related orbits that violate this assumption. We were able to test how 
well the assumption holds for our mock halo stars. Perhaps surpris¬ 
ingly, we do not find any systematic correlation of the recovered 
halo mass with the infall redshift of tracer subsamples. This sug¬ 
gests that the dynamical state of halo tracers depends on other fac¬ 
tors, such as their orbits, and not only their infall time. Dynamical 
relaxation is nevertheless a factor: excluding stars stripped from 
surviving satellites improves the agreement between best-fit and 
true halo masses in two cases (halos D and E). This cut eliminates 
dynamically hot structures that can be identified by eye in these 
halos. 

Beyond all these assumptions and uncertainties in the model 
itself, in real observations the maximum observable radius of dy¬ 
namical tracers may be much smaller than the halo virial radius. We 
found tracer subsamples selected over different ranges of radius can 
give significantly different estimates of the host halo mass, even 
if the three parameters describing the density of tracers are fixed 
to be those derived from the whole tracer population. An outer 
radius limit results in biased measurements of M 200 if it signif¬ 
icantly smaller than the virial radius. For example, the recovered 
halo masses of halos A, B, C and D converge for outer radius limits 
larger than r ~ 0.4f?2oo but give systematically larger masses for 
smaller radial limits. For one halo, E. this overestimation occurs 
for limits r < 0.8f?200- There are two reasons behind this radial 
dependence: stars at different radii have significantly different dy¬ 
namical state and extrapolations to larger and smaller radii become 
less accurate when only a limited radial range is sampled. 

Real surveys have complex selection functions for stars, which 
depend on both radial distance and angular position. Particular 
classes of tracers may be very sparsely sampled. The observed par¬ 
allax, radial and tangential velocities of halo stars include observa¬ 
tional uncertainties which depend strongly on distance. Although 
a large sample of tracers with exact coordinate and velocity infor¬ 
mation from our mock stellar halo catalogue have enabled us to 
investigate the model performance, it will be important to consider 
realistic observational errors and sample selection effects in future 
studies aimed at forecasting the performance of real surveys. 

We conclude that methods to estimate the mass of the Milky 
Way halo using the kind of distribution functions we have investi¬ 
gated here need to be used with extreme caution. This is particularly 
true when estimating the total virial mass. Restricting the estimate 
to the mass interior to r ~ 0 . 2 f ?200 is considerably more reliable. 
In any case, mock catalogues like th ose w e hav e analysed here and 
made publicly available in lLowing et al.1 ( 12015 4 are required to as¬ 
sess the reliability of any particular mass estimation method. 
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APPENDIX A: ORIGIN OF THE BIAS IN /3 
Al Origin of the bias 

In Fig.[ 6 ]we showed that there is a systematic bias between the best- 
fit and true value of /3. We now show in the left panels of Fig. IA1I 
the phase-space distribution of stars in halo A (top panel, binding 
energy, E, versus angular momentum, L) and the one dimensional 
angular momentum distribution at two fixed values of E (the sec¬ 
ond and third panels from the top, respectively) indicated by the 
horizontal dashed lines in the top paneljj E and L have been nor¬ 
malised so that they are dimensionless (see Eqn. [13. We only show 
results based on halo A; for the other halos our conclusions are the 
same. 

In the top panel we see that, for a fixed value of E, there is 
an upper limit to L which increases with decreasing E. This is the 
maximum allowed value of angular momentum, corresponding to 
circular orbits with zero radial velocity at fixed E. In the next two 
panels, we show the angular momentum distributions at different 
values of E have similar features. They are both flat at small val¬ 
ues of L and drop quickly when L approaches its upper limit. For 
comparison, we also plot two lines of the form F(E, L) oc L -2/3 , 
where /3 is the velocity anisotropy obtained from stars in the en¬ 
ergy slice, /3b (red dashed lines), or the full sample, /3 OV er a ii (green 
dashed lines). /3 0 veraii is the same in both panels and is simply the 
true value of /3 from Table |T| (0.66). The values of /3e are 0.442 
and 0.618 for the middle and bottom panels respectively. Neither 
of the two could give a satisfactory description of the true distri¬ 
bution (blue curve). This implies that the physical interpretation of 
the power-law index in our distribution function as —2/3 is inac¬ 
curate, and this is the origin for the systematic bias of /3. The true 
distribution function must be more complex. 

We also notice that the best-fit value of /3 for halo A is 0.458 
(Table QJ, which predicts a power-law slope shallower than the 
green dashed lines and is also shallower than the red dashed line 
in the bottom panel, but it is still a poor match to the distributions 
(blue curve) in Fig. lAll If we fix the power-law slope in the model 
according to the true anisotropy of the full sample, this results in 
better agreement with tangential velocity distribution but a much 
poorer agreement with the radial velocity distribution. Afterall, our 
maximum likelihood approach is designed to fit the velocity and 
spatial distributions of stars, not the distributions of binding energy 
or angular momentum. 


' Note the quantity we are plotting here is F{L) = F(L\E). To obtain 
this distribution empirically one has to properly account for the density of 
state in (E,L) space (see|Woitak et al. 120081 for more details). 


A2 Why stars are more radially biased than dark matter? 

The /3 profile of dark matter particles have bee n studied in earlier 
works. For example, IWoitak et al] d2008l , l2009t) looked at the dis¬ 
tribution functions of dark matter particles in halos of mass 10 14 
to 10 15 Mq. Although the details of their modelling and the mass 
range of halos are different from ours, their model distribution func¬ 
tion can recover well the true /3 of dark matter particles in their sim¬ 
ulation. We therefore examine the angular momentum distribution 
of dark matter particles in our simulations in the three right panels 
of Fig.lAll 

In the top panel, we see dark matter particles can extend to 
much lower binding energy than stars. Black curves in the mid¬ 
dle and bottom panels show the angular momentum distribution 
at two fixed values of E. We again plot two lines of the form 
F(E, L) oc L~ 213 . Red and green lines are predicted from the ve¬ 
locity anisotropy of dark matter particles in the energy slice or from 
the the full sample respectively. At E ~ 10 -o ' 4 u 2 , the agreement 
between the red dashed lines and the shape of the L distributions is 
quite poor, whereas at a lower binding energy (E ~ 10 _11 i;J), we 
see a better agreement. We have looked at many different choices of 
E in this regard, and found that for less bound dark matter particles, 
their velocity anisotropy correctly predicts the power-law slope of 
their L distribution. This means the model distribution function de¬ 
scribes better systems of less bound dark matter particles. However, 
for dark matter particles that are more tightly bound, the velocity 
anisotropy is not as well correlated with the power-law slope of 
the L distribution. This is the same as the stellar case, although the 
discrepancy for dark matter particles is smaller. 

Stars in the stellar halo are clearly a biased population of trac¬ 
ers with respect to dark matter particles in the simulation. Their 
orbits are more radial (Fig. [3 corresponding to a higher /3. How¬ 
ever, the difference in /3 is not only because stars are more dynam¬ 
ically bound than dark matter particles: we have explicitly checked 
that, for a given fraction of the most bound dark matter particles, 
orbits are still more tangentially biased than stars with the same 
range of binding energy. The fact that stars are more radially bi¬ 
ased than dark matter particles thus has more fundamental physical 
origin. First of all, in our model halo stars are all accreted from 
subhalos, while dark matter particles are added to the main halo by 
both clumpy and smooth accretion. We have calculated the velocity 
anisotropy of dark matter accreted from subhalos only, and found 
these particles are more radially biased than all the dark matter par¬ 
ticles as a whole. This is probably because the clumps in which 
these particles are accreted (i.e. subhalos) have more radially biased 
orbits. Furthermore, halo stars in our analysis are tags placed on 
the most bound dark matter particles in progenitor su bhalo s, which 
have then been stripped and mixed into the main halo. [Lowing et al] 
d2015t) have found the halo stars are dominated by contributions 
from a few massive satellites. As the most bound parts of these 
satellites have been stripped into the halo, the satellites are more 
likely to have been on highly radial orbits, imparting a radial bias 
to halo stars. In contrast, dark matter particles enter the main halos 
in our simulations through quite different mechanisms, with both 
clumpy and smooth accretion dWang et~al]|201 ll) . 
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Figure Al. phase-space distributions of mock stars in halo A (left) and of dark matter in the same halo (right). The top panels in either plots show 2D 
distributions in the plane of binding energy, E, and angular momentum, L. E and L have been scaled to be dimensionless. The middle and bottom panels 
show ID angular momentum distributions of stars (likewise DM) at two fixed values of E, indicated by the horizontal black dashed curves in the top panels. 
Red and green dashed lines are power-law distributions with arbitrary normalizations, predicted from the velocity anisotropy of all stars (or DM) and stars (or 
DM) in each binding energy range. 


APPENDIX B: UNCERTAINTIES IN MODELLING THE 
HALO BOUNDARIES 

For all the analysis in the main text, we have been assuming the spa¬ 
tial extent of both NFW halos and tracers are infinite (r max ,h = oo 
and r max ,t = oo). It is, however, necessary to investigate whether 
the different choices of halo and tracer boundaries could affect our 
measured halo properties. 

As we have mentioned in Sec. 13.21 in principle tracer bound¬ 
aries (r maX! t) can be different from halo boundaries. Here for sim¬ 
plicity we assume r maXi h = r max , t . We tried four different choices 
of halo boundaries (r max ,h), ranging from twice to five times the 
halo virial radius. We avoid using boundaries at exactly the halo 
virial radius because our mock halo stars can extend beyond R 200 , 
while the mass distribution in the FoF group distribute continuously 
and extend further than R 200 ■ A sharp cut at R 200 is thus not real¬ 
istic. 

The best-fit host halo masses and concentrations as functions 
of halo boundaries are presented in Fig. lBlI The velocity anisotropy 
/3 almost does not change with the different choices of halo bound¬ 
aries, and thus we do not show them. Dashed red lines are true 
values of halo masses and concentrations. 

The measured halo masses increase with the decrease in halo 
boundaries, and the halo concentrations decrease accordingly, re¬ 
flecting again the strong correlation between the two parameters. 
The choice of halo boundary that gives the best match between 
measured and true halo mass varies from halo to halo. For halo B, 
the best fit halo masses and concentrations almost do not change 
with the choice of halo boundaries when r max ,h ^ 3f?200 and 


agree well with the true values. At r maXj h = 2 R 200 , the measured 
host halo mass gets significantly larger, indicating for halo B fi¬ 
nite halo boundaries do not help to improve the fitting. For halo 
C, r maXj h = 2 R 200 gives a very good match between the best-fit 
and true halo mass and concentration, demonstrating a finite halo 
boundary works better than infinite boundaries at least for halo C. 

The estimated halo mass of halo A is closest to the true value 
when halo boundary is chosen at twice the virial radius. However, 
the estimated halo concentration at that boundary deviates signifi¬ 
cantly from the true concentration, suggesting the discrepancy be¬ 
tween best-fit and true host halo mass of halo A could not be domi¬ 
nated by how boundaries are modeled. For halo D and E, we know 
already because of the existence of unrelaxed dynamical structures, 
the host halo masses are significantly underestimated, and thus the 
best match between measured and true halo parameters at twice 
(halo E) and four times (halo D) the virial radius demonstrates the 
entangling of different model uncertainties, which canceled with 
each other to give a good prediction. 
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Figure Bl. The measured host halo mass (top row) and concentration (bottom row) as functions of the halo boundaries in the model. 


APPENDIX C: DARK MATTER PARTICLES AS 
TRACERS 

In Sec.[A]we show the velocity anisotropies of dark matter particles 
agree better with the distribution function model than those of stars. 
Furthermore, dark matter particles are more radially extended than 
stars and might probe better the underlying potential in outskirts. 
Thus we ask whether better constraints on the halo properties can 
be achieved by using dark matter particles as tracers. Obviously, 
it is not possible to directly observe the dynamics of dark matter, 
but asking this question helps to deepen our understanding of the 
model. The answer is, unfortunately, no. By using dark matter par¬ 
ticles as tracers we end up with significant overestimates of the 
host halo mass, at least for our five Milky Way analog halos. These 
measurements are shown in Table ed where we have used a ran¬ 
domly selected subsample of all dark matter particles in the halo 
FoF group (one particle out of every 5,000 in the simulation). We 
have explicitly checked that this conclusion does not change if we 
randomly select different subsamples of dark matter, remove dark 
matter particles in substructures or restrict them to be inside the 
halo virial radius. Both radial and tangential velocities have been 
used in this analysis. 

To explore the reasons behind this, we present in Fig. eh 
phase-space contour plots for dark matter particles and stars in the 
simulation and compare these with realizations drawn directly from 
the model distribution function. We only show plots based on halo 
A; for the other halos the conclusion is the same. Distributions of 
binding energy versus radial velocity, v r , tangential velocity, vt, 
and radius, r, have been plotted separately, so that we are able to 
see how well the model prediction agrees with the true distribution 
of v r , vt and r for dark matter particles in the simulation. 

It is very clear to see that, with true halo parameters, the model 
predictions deviate significantly from the empirical distribution of 
r, v r and vt at low binding energy (left column). On the other hand. 


the best fit model agrees much better with the data (middle col¬ 
umn). This improved agreement is caused by an overestimate of the 
host halo mass, leading to a deeper potential and increased binding 
energy for tracer particles. As a result, the sample becomes more 
dynamically bound and agrees better with the model. 

Our conclusion is therefore that, although the model distribu¬ 
tion function gives a better approximation in the velocity anisotropy 
of dark matter particles, the predicated phase-space distribution at 
the low binding energy end is very poor. By construction, the dis¬ 
tribution function only describes closed systems. This not only re¬ 
quires that the tracer population should be bound and truncated, but 
also results in a sharp cut-off in binding energy that is a poor de¬ 
scription of particle distribution in the low energy end. Compared 
with all dark matter particles, the star particles in our simulation are 
more dynamically bound and centrally concentrated (see the right 
column) and are thus not sensitive to the low binding energy tail of 
the dark matter distribution. Using stars rather than dark matter par¬ 
ticles therefore insulates the fit from deficiencies of the distribution 
function model at the low binding energy tail. 
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Table Cl. best-fit parameters from the phase-space distribution of dark matter particles 



A 

B 

C 

D 

E 

M 200 [io 12 M Q ] 

C 200 
r s [kpc] 

logiQ p a [M 0 /kpc 3 ] 
R200 [kpc] 

P 

2.811 ±0.024 
4.458 ± 0.098 
63.499 ± 1.116 
5.997 ± 0.021 
283.078 ± 7.965 
0.266 ± 0.004 

1.227 ± 0.012 
4.845 ±0.105 
44.315 ±0.773 
6.078 ± 0.021 
214.715 ± 5.956 
0.159 ±0.005 

2.868 ± 0.024 

6.031 ±0.106 
47.257 ± 0.683 
6.297 ±0.018 
285.003 ± 6.482 
0.259 ± 0.005 

2.727 ± 0.029 
3.903 ±0.104 
71.799 ± 1.516 
5.868 ± 0.025 
280.225 ± 9.529 
0.217 ± 0.006 

1.776 ±0.015 
5.296 ±0.108 
45.867 ± 0.757 
6.166 ±0.020 
242.901 ± 6.357 
0.090 ± 0.005 


DM:true DM: best fit star: true 



I£[(km/s) 2 ] ^[(km/s) 2 ] i£[(km/s) 2 ] 


Figure Cl. The phase-space density of dark matter particles or stars in halo 
A (green solid) and model predictions (red dashed). The contours mark the 
10th, 30th, 60th and 90th percentiles of the 2D density distribution in pa¬ 
rameter plane. We present contour plots of binding energy, E, versus radius, 
r, radial velocity, v r , and tangential velocity, vt. For simplicity we only use 
the magnitudes of v r and vt , so all quantities are positive. In deducing the 
binding energy, we use the analytical NFW potential model. Green con¬ 
tours in the left and middle columns are based on dark matter particles in 
the simulation, while in the right column we plot contours for stars. For the 
left, middle and right columns, true halo parameters, dynamical best-fit halo 
parameters from dark matter particles and true halo parameters are adopted 
in the potential model respectively. The lines are iso-density contours that 
contain the 10,30,60 and 90% densest cells. 
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